Marrow FACS Notebook

Specify the tissue of interest, run the boilerplate code which sets up the functions and environment, load the tissue object.

tissue_of_interest = "Marrow"
library(here)

here() starts at /Users/olgabot/code/tabula-muris

source(here("00_data_ingest", "02_tissue_analysis_rmd", "boilerplate.R"))

Loading required package: ggplot2
Loading required package: cowplot

Attaching package: ‘cowplot’

The following object is masked from ‘package:ggplot2’:

    ggsave

Loading required package: Matrix

tiss = load_tissue_facs(tissue_of_interest)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

[1] "Scaling data matrix"

  |                                                                                                  
  |                                                                                            |   0%
  |                                                                                                  
  |====                                                                                        |   4%
  |                                                                                                  
  |========                                                                                    |   8%
  |                                                                                                  
  |============                                                                                |  12%
  |                                                                                                  
  |===============                                                                             |  17%
  |                                                                                                  
  |===================                                                                         |  21%
  |                                                                                                  
  |=======================                                                                     |  25%
  |                                                                                                  
  |===========================                                                                 |  29%
  |                                                                                                  
  |===============================                                                             |  33%
  |                                                                                                  
  |==================================                                                          |  38%
  |                                                                                                  
  |======================================                                                      |  42%
  |                                                                                                  
  |==========================================                                                  |  46%
  |                                                                                                  
  |==============================================                                              |  50%
  |                                                                                                  
  |==================================================                                          |  54%
  |                                                                                                  
  |======================================================                                      |  58%
  |                                                                                                  
  |==========================================================                                  |  62%
  |                                                                                                  
  |=============================================================                               |  67%
  |                                                                                                  
  |=================================================================                           |  71%
  |                                                                                                  
  |=====================================================================                       |  75%
  |                                                                                                  
  |=========================================================================                   |  79%
  |                                                                                                  
  |=============================================================================               |  83%
  |                                                                                                  
  |================================================================================            |  88%
  |                                                                                                  
  |====================================================================================        |  92%
  |                                                                                                  
  |========================================================================================    |  96%
  |                                                                                                  
  |============================================================================================| 100%

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Visualize top genes in principal components

Later on (in FindClusters and TSNE) you will pick a number of principal components to use. This has the effect of keeping the major directions of variation in the data and, ideally, supressing noise. There is no correct answer to the number to use, but a decent rule of thumb is to go until the plot plateaus.

PCElbowPlot(object = tiss)

Choose the number of principal components to use.

# Set number of principal components. 
n.pcs = 18

The clustering is performed based on a nearest neighbors graph. Cells that have similar expression will be joined together. The Louvain algorithm looks for groups of cells with high modularity–more connections within the group than between groups. The resolution parameter determines the scale…higher resolution will give more clusters, lower resolution will give fewer.

For the top-level clustering, aim to under-cluster instead of over-cluster. It will be easy to subset groups and further analyze them below.

# Set resolution 
res.used <- 0.5
tiss <- FindClusters(object = tiss, reduction.type = "pca", dims.use = 1:n.pcs, 
    resolution = res.used, print.output = 0, save.SNN = TRUE)

To visualize

# If cells are too spread out, you can raise the perplexity. If you have few cells, try a lower perplexity (but never less than 10).
tiss <- RunTSNE(object = tiss, dims.use = 1:n.pcs, seed.use = 10, perplexity=84)

# note that you can set do.label=T to help label individual clusters
TSNEPlot(object = tiss, do.label = T)

Compare to previous annotations

method = "facs"
tiss = compare_previous_annotation(tiss, tissue_of_interest, method, filename = here('00_data_ingest', '03_tissue_annotation_csv', 
                    paste0(tissue_of_interest, "_", method, "_annotation_rahul.csv")))

                        Cd3e+_Klrb1+ B cell                  common lymphoid progenitor 
                                         37                                         129 
                      Dntt- late pro-B cell                       Dntt+ late pro-B cell 
                                        203                                          92 
       granulocyte monocyte progenitor cell   hematopoietic stem cell enriched fraction 
                                         65                                          50 
                      HSC enriched fraction                immature natural killer cell 
                                        158                                          31 
                         immature NK T cell                             immature T cell 
                                         46                                          53 
                                  MAA001844                  mature natural killer cell 
                                        194                                          54 
    megakaryocyte-erythroid progenitor cell                                          NA 
                                         46                                         140 
                    pre-natural killer cell                           regulatory T cell 
                                         19                                          32 
Slamf1-negative multipotent progenitor cell Slamf1-positive multipotent progenitor cell 
                                        548                                          89 
                                             
                                                0   1   2   3   4   5   6   7   8   9  10  11  12
  Cd3e+_Klrb1+ B cell                           0   0   0   0   0   0   0   0   0   0   0  33   4
  common lymphoid progenitor                  127   0   0   0   2   0   0   0   0   0   0   0   0
  Dntt- late pro-B cell                         0   0   0   0   0   0   0   0 203   0   0   0   0
  Dntt+ late pro-B cell                         0   0   0   0   0   0   0   0  92   0   0   0   0
  granulocyte monocyte progenitor cell         62   0   0   0   3   0   0   0   0   0   0   0   0
  hematopoietic stem cell enriched fraction    45   0   0   0   5   0   0   0   0   0   0   0   0
  HSC enriched fraction                        93   0   0   0  65   0   0   0   0   0   0   0   0
  immature natural killer cell                  0   0   0   0   0   0   0   0   0   0   0   2  29
  immature NK T cell                            0   0   0   0   0   0   0   0   0   0   0  45   1
  immature T cell                               0   0   0   0   0   0   0   0   0   0   0  53   0
  MAA001844                                     0   0   0   0 194   0   0   0   0   0   0   0   0
  mature natural killer cell                    0   0   0   0   0   0   0   0   0   0   0   2  52
  megakaryocyte-erythroid progenitor cell      46   0   0   0   0   0   0   0   0   0   0   0   0
  NA                                          106   0   0   0   0   2   0  11   6  11   0   2   2
  pre-natural killer cell                       0   0   0   0   0   0   0   0   0   0   0   3  16
  regulatory T cell                             0   0   0   0   0   0   0   0   0   0   0  32   0
  Slamf1-negative multipotent progenitor cell 446   0   0   0 102   0   0   0   0   0   0   0   0
  Slamf1-positive multipotent progenitor cell  89   0   0   0   0   0   0   0   0   0   0   0   0

                            basophil                          granulocyte 
                                  24                                  758 
granulocyte monocyte progenitor cell              granulocytopoietic cell 
                                 134                                  211 
                     immature B cell                           macrophage 
                                 365                                  174 
                            monocyte                                   NA 
                                 257                                  140 
                        naive B cell                           pre-B cell 
                                 676                                  517 
                                      
                                         0   1   2   3   4   5   6   7   8   9  10  11  12
  basophil                               0   0   0   0   0   0   0   0   0   0   0   0  24
  granulocyte                            0   0   0 390   0 368   0   0   0   0   0   0   0
  granulocyte monocyte progenitor cell  62   0   0   0   3   0   0  69   0   0   0   0   0
  granulocytopoietic cell                0   0   0   0   0   1   0   0   0 210   0   0   0
  immature B cell                        0  21   3   0   0   0 341   0   0   0   0   0   0
  macrophage                             0   0   0   0   1   0   0   0   0   0 173   0   0
  monocyte                               1   0   0   0   0   0   0 256   0   0   0   0   0
  NA                                   106   0   0   0   0   2   0  11   6  11   0   2   2
  naive B cell                           0 676   0   0   0   0   0   0   0   0   0   0   0
  pre-B cell                             0   0 514   0   0   0   3   0   0   0   0   0   0

TSNEPlot(object = tiss, do.return = TRUE, group.by = "previous_cell_ontology_class")

TSNEPlot(object = tiss, do.return = TRUE, group.by = "previous_free_annotation")

Check expression of genes of interset.

Dotplots let you see the intensity of exppression and the fraction of cells expressing for each of your genes of interest.

How big are the clusters?

table(tiss@ident)

   0    1    2    3    4    5    6    7    8    9   10   11   12 
1015  697  517  390  372  371  344  336  301  221  173  172  128 

Assigning cell type identity to clusters

At a coarse level, we can use canonical markers to match the unbiased clustering to known cell types:

previous_annotation_table(tiss)

[1] "previous_free_annotation"
                                             
                                                0   1   2   3   4   5   6   7   8   9  10  11  12
  Cd3e+_Klrb1+ B cell                           0   0   0   0   0   0   0   0   0   0   0  33   4
  common lymphoid progenitor                  127   0   0   0   2   0   0   0   0   0   0   0   0
  Dntt- late pro-B cell                         0   0   0   0   0   0   0   0 203   0   0   0   0
  Dntt+ late pro-B cell                         0   0   0   0   0   0   0   0  92   0   0   0   0
  granulocyte monocyte progenitor cell         62   0   0   0   3   0   0   0   0   0   0   0   0
  hematopoietic stem cell enriched fraction    45   0   0   0   5   0   0   0   0   0   0   0   0
  HSC enriched fraction                        93   0   0   0  65   0   0   0   0   0   0   0   0
  immature natural killer cell                  0   0   0   0   0   0   0   0   0   0   0   2  29
  immature NK T cell                            0   0   0   0   0   0   0   0   0   0   0  45   1
  immature T cell                               0   0   0   0   0   0   0   0   0   0   0  53   0
  MAA001844                                     0   0   0   0 194   0   0   0   0   0   0   0   0
  mature natural killer cell                    0   0   0   0   0   0   0   0   0   0   0   2  52
  megakaryocyte-erythroid progenitor cell      46   0   0   0   0   0   0   0   0   0   0   0   0
  NA                                          106   0   0   0   0   2   0  11   6  11   0   2   2
  pre-natural killer cell                       0   0   0   0   0   0   0   0   0   0   0   3  16
  regulatory T cell                             0   0   0   0   0   0   0   0   0   0   0  32   0
  Slamf1-negative multipotent progenitor cell 446   0   0   0 102   0   0   0   0   0   0   0   0
  Slamf1-positive multipotent progenitor cell  89   0   0   0   0   0   0   0   0   0   0   0   0
[1] "previous_cell_ontology_class"
                                      
                                         0   1   2   3   4   5   6   7   8   9  10  11  12
  basophil                               0   0   0   0   0   0   0   0   0   0   0   0  24
  granulocyte                            0   0   0 390   0 368   0   0   0   0   0   0   0
  granulocyte monocyte progenitor cell  62   0   0   0   3   0   0  69   0   0   0   0   0
  granulocytopoietic cell                0   0   0   0   0   1   0   0   0 210   0   0   0
  immature B cell                        0  21   3   0   0   0 341   0   0   0   0   0   0
  macrophage                             0   0   0   0   1   0   0   0   0   0 173   0   0
  monocyte                               1   0   0   0   0   0   0 256   0   0   0   0   0
  NA                                   106   0   0   0   0   2   0  11   6  11   0   2   2
  naive B cell                           0 676   0   0   0   0   0   0   0   0   0   0   0
  pre-B cell                             0   0 514   0   0   0   3   0   0   0   0   0   0

# stash current cluster IDs
tiss <- StashIdent(object = tiss, save.name = "cluster.ids")
# enumerate current cluster IDs and the labels for them
cluster.ids <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)
free_annotation <- c(NA, NA, "pre-B cell (Philadelphia nomenclature)",
                     NA, NA, NA, NA, NA, NA, NA, NA, NA, NA)
cell_ontology_class <-c(
  NA, # Cluster 0 is a mix of progenitors and precursors
  "naive B cell", 
  "precursor B cell", 
  "granulocyte", 
  "hematopoietic precursor cell", 
  "granulocyte",
  "immature B cell", 
  NA, # Cluster 7 is a mix of monocytes and granulocyte progenitors
  "late pro-B cell", 
  "granulocytopoietic cell",
  "macrophage",
  NA, # Cluster 11 is a mix of T, NKT and B cells
  NA # Cluster 12 is a mix of basophil, B cells, immature NK cells, mature NK cells, pre-natural killer cells
)
tiss = stash_annotations(tiss, cluster.ids, free_annotation, cell_ontology_class)
data.frame(cluster.ids, cell_ontology_class, free_annotation)

TSNEPlot(object = tiss, do.label = TRUE, pt.size = 0.5, group.by='cell_ontology_class')

Checking for batch effects

Color by metadata, like plate barcode, to check for batch effects.

TSNEPlot(object = tiss, do.return = TRUE, group.by = "plate.barcode")

Print a table showing the count of cells in each identity category from each plate.

table(as.character(tiss@ident), as.character(tiss@meta.data$plate.barcode))

    
     D041912 D042044 D042193 D042467 D042479 D045139 D045140 MAA000409 MAA000596 MAA000600 MAA000639
  0        1       0       2       1       0       0       0         0       290       333         0
  1        2      89      96      39      53       0       0         2         1         5         0
  10       2      31       2       1      15       0       0         0         0         0         0
  11       0       0      50      20       2       0       0         0         0         0         0
  12       6       0      33      19       1       3       2         1         0         0         0
  2        0      53      45      16      59       0       0         0         0         0         0
  3       68       2       3       2       1      29      84       132         0         0        62
  4        0       0       0       0       0       0       2         0         0         1         0
  5       75       0       2       1       0      38      59        75         0         0       116
  6        0      42      48      18      33       0       0         0         1         0         0
  7       49       1       8       4       3      45      68        50         0         0        45
  8        0      89       7       3      82       0       0         0         0         0         0
  9       69       0       2       0       0      14      16        62         0         0        50
    
     MAA000652 MAA000844 MAA000848 MAA001842 MAA001844 MAA001883 MAA001884 MAA001887 MAA001888
  0          0         1       346         0         3         0         0        38         0
  1         52        67         2        38         7        38       171         5        30
  10        38         9         0         1         0        26        37         0        11
  11         0        47         0         1         2         1         2         1        46
  12         0        23         0         3         1         1         0         2        33
  2        122        62         0         4         0       104        49         0         3
  3          2         5         0         0         0         0         0         0         0
  4          0         0         3         0       190         0         1       174         1
  5          2         2         0         0         0         0         0         0         1
  6         36        42         0         5         2        45        64         3         5
  7          0         8         0         2        28         0         0        17         8
  8         46         5         0         1         1        33        26         7         1
  9          1         2         0         0         3         0         0         2         0

Redo Subcluster as per the code from pancreas plate Rmd, which is without regressing out ribosomal genes

Subset A == T_NK_NKT cell (Clusters 11, 12)

subtissA = SubsetData(tiss, ident.use = c(11, 12))

subtissA <- subtissA %>% ScaleData() %>% 
  FindVariableGenes(do.plot = TRUE, x.high.cutoff = Inf, y.cutoff = 0.5) %>%
  RunPCA(do.print = FALSE)

[1] "Scaling data matrix"

  |                                                                                                  
  |                                                                                            |   0%
  |                                                                                                  
  |============================================================================================| 100%

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

PCHeatmap(object = subtissA, pc.use = 1:3, cells.use = 100, do.balanced = TRUE, label.columns = FALSE, num.genes = 8)

PCElbowPlot(subtissA)

sub.n.pcs = 4
sub.res.use = 1
subtissA <- subtissA %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, 
    resolution = sub.res.use, print.output = 0, save.SNN = TRUE) %>%
    RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=44)
TSNEPlot(object = subtissA, do.label = T, pt.size = 1.2, label.size = 4)

Check expression of genes of interset.

{r} # subtissA.markers <- FindAllMarkers(object = subtissA, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) #

{r} # subtissA.markers %>% group_by(cluster) %>% top_n(6, avg_logFC) #

TSNEPlot(object = subtissA, do.return = TRUE, group.by = "previous_cell_ontology_class")

TSNEPlot(object = subtissA, do.return = TRUE, group.by = "previous_free_annotation")

previous_annotation_table(subtissA)

[1] "previous_free_annotation"
                              
                                0  1  2  3  4  5  6
  Cd3e+_Klrb1+ B cell           0  0 37  0  0  0  0
  immature natural killer cell  0  0  0  0 31  0  0
  immature NK T cell            1  0  6 36  3  0  0
  immature T cell              50  0  0  3  0  0  0
  mature natural killer cell    0 48  2  0  3  1  0
  NA                            0  0  0  1  0  0  3
  pre-natural killer cell       0  0  0  0  0  0 19
  regulatory T cell            31  0  0  1  0  0  0
[1] "previous_cell_ontology_class"
          
            0  1  2  3  4  5  6
  basophil  0  0  0  0  0 24  0
  NA        0  0  0  1  0  0  3

The multitude of clusters of each type correspond mostly to individual animals/sexes.

table(FetchData(subtissA, c('mouse.id','ident')) %>% droplevels())

        ident
mouse.id  0  1  2  3  4  5  6
  3_10_M 27  8  7 10  9  4  6
  3_38_F  1  0  3  1  5  1  0
  3_39_F 12  8 13 21 15 11  7
  3_8_M  30 18 16  8  7  5  8
  3_9_M  12 14  6  1  1  4  1

sub.cluster.ids <- c(0, 1, 2, 3, 4, 5, 6)
sub.free_annotation <- c("regulatory and immature T cell", NA, "Cd3e+ Klrb1+ B cell",  NA, NA, NA, NA)
""

[1] ""

sub.cell_ontology_class <-c("T cell", "mature natural killer cell","B cell",  "immature NK T cell", "immature natural killer cell", "basophil",  "pre-natural killer cell")
subtissA = stash_annotations(subtissA, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class)
tiss = stash_subtiss_in_tiss(tiss, subtissA)
head(tiss@meta.data)

data.frame(sub.cluster.ids, sub.cell_ontology_class, sub.free_annotation)

Checking for batch effects

Color by metadata, like plate barcode, to check for batch effects.

TSNEPlot(object = subtissA, do.return = TRUE, group.by = "plate.barcode")

Final coloring

Color by cell ontology class on the original TSNE.

TSNEPlot(object = subtissA, do.label = TRUE, do.return = TRUE, group.by = "cell_ontology_class")

TSNEPlot(object = subtissA, do.label = TRUE, do.return = TRUE, group.by = "free_annotation")

Subset B == regulatory_immature T cell (Cluster 0)

Multiple cell-types in the cluster 0 of subtiss (original cluster 7)

subcluster the cluster 0 (regulatory_immature T cell) of subtiss (original cluster 7)

subtissB = SubsetData(subtissA, ident.use = c(0))

subtissB <- subtissB %>% ScaleData() %>% 
  FindVariableGenes(do.plot = TRUE, x.high.cutoff = Inf, y.cutoff = 0.5) %>%
  RunPCA(do.print = FALSE)

[1] "Scaling data matrix"

  |                                                                                                  
  |                                                                                            |   0%
  |                                                                                                  
  |============================================================================================| 100%

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

PCHeatmap(object = subtissB, pc.use = 1:3, cells.use = 150, do.balanced = TRUE, label.columns = FALSE, num.genes = 8)

PCElbowPlot(subtissB)

sub.n.pcs = 3
sub.res.use = 1
subtissB <- subtissB %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, 
    resolution = sub.res.use, print.output = 0, save.SNN = TRUE) %>%
    RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=12)
TSNEPlot(object = subtissB, do.label = T, pt.size = 1.2, label.size = 4)

previous_annotation_table(subtissB)

[1] "previous_free_annotation"
                    
                      0  1  2
  immature NK T cell  0  0  1
  immature T cell    39  1 10
  regulatory T cell   1 25  5
[1] "previous_cell_ontology_class"
< table of extent 0 x 3 >

Check expression of genes of interset.

{r} # subtissB.markers <- FindAllMarkers(object = subtissB, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) #

{r} # subtissB.markers %>% group_by(cluster) %>% top_n(6, avg_logFC) #

From these genes, it appears that the clusters represent:

The multitude of clusters of each type correspond mostly to individual animals/sexes.

table(FetchData(subtissB, c('mouse.id','ident')) %>% droplevels())

        ident
mouse.id  0  1  2
  3_10_M 10  9  8
  3_38_F  1  0  0
  3_39_F  0  9  3
  3_8_M  23  5  2
  3_9_M   6  3  3

previous_annotation_table(subtissB)

[1] "previous_free_annotation"
                    
                      0  1  2
  immature NK T cell  0  0  1
  immature T cell    39  1 10
  regulatory T cell   1 25  5
[1] "previous_cell_ontology_class"
< table of extent 0 x 3 >

sub.cluster.ids <- c(0, 1, 2)
sub.free_annotation <-c(NA, NA, NA)
sub.cell_ontology_class <-c("immature T cell", "regulatory T cell", "immature T cell")
subtissB = stash_annotations(subtissB, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class)
tiss = stash_subtiss_in_tiss(tiss, subtissB)
data.frame(sub.cluster.ids, sub.cell_ontology_class, sub.free_annotation)

Checking for batch effects

Color by metadata, like plate barcode, to check for batch effects.

TSNEPlot(object = subtissB, do.return = TRUE, group.by = "plate.barcode")

Final coloring

Color by cell ontology class on the original TSNE.

TSNEPlot(object = subtissB, do.label = TRUE, do.return = TRUE, group.by = "cell_ontology_class")

Subset C == late pro-B cell (Cluster 8)

late pro-B cells are reportedly committed to B cell lineage as they begin expressing Pax5. They also have Rag1, Rag2, and Dntt expression. At single cell resolution there seem to be two subsets of late pro-B cells: Dntt+ and Dntt-.

It appears there are multiple novel cell-types within the Pax5+ late pro-B cell

subtissC = SubsetData(tiss, ident.use = c(8))

subtissC <- subtissC %>% ScaleData() %>% 
  FindVariableGenes(do.plot = TRUE, x.high.cutoff = Inf, y.cutoff = 0.5) %>%
  RunPCA(do.print = FALSE)

[1] "Scaling data matrix"

  |                                                                                                  
  |                                                                                            |   0%
  |                                                                                                  
  |============================================================================================| 100%

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

PCHeatmap(object = subtissC, pc.use = 1:3, cells.use = 500, do.balanced = TRUE, label.columns = FALSE, num.genes = 8)

PCElbowPlot(subtissC)

sub.n.pcs = 4
sub.res.use = 1
subtissC <- subtissC %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, 
    resolution = sub.res.use, print.output = 0, save.SNN = TRUE) %>%
    RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=42)
TSNEPlot(object = subtissC, do.label = T, pt.size = 1.2, label.size = 4)

previous_annotation_table(subtissC)

[1] "previous_free_annotation"
                       
                         0  1  2  3  4
  Dntt- late pro-B cell 85 55 63  0  0
  Dntt+ late pro-B cell  1 15  0 44 32
  NA                     0  2  1  3  0
[1] "previous_cell_ontology_class"
    
     0 1 2 3 4
  NA 0 2 1 3 0

Check expression of genes of interset.

{r} # subtissC.markers <- FindAllMarkers(object = subtissC, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) #

{r} # subtissC.markers %>% group_by(cluster) %>% top_n(6, avg_logFC) #

genes_to_check = c('Dntt')
# FeaturePlot(subtissC, genes_to_check, pt.size = 1, nCol = 3)
DotPlot(subtissC, genes_to_check, col.max = 2.5, plot.legend = T, do.return = T) + coord_flip()

The multitude of clusters of each type correspond mostly to individual animals/sexes.

table(FetchData(subtissC, c('mouse.id','ident')) %>% droplevels())

        ident
mouse.id  0  1  2  3  4
  3_10_M 24  3 17  4  3
  3_38_F  9  2  7  2  8
  3_39_F 11  9  4 14  3
  3_8_M  16 25 22 19  7
  3_9_M  26 33 14  8 11

sub.cluster.ids <- c(0, 1, 2, 3, 4)
sub.free_annotation <-
  c(
  "Dntt- late pro-B cell",
  "Dntt- late pro-B cell",
  "Dntt- late pro-B cell",
  "Dntt+ late pro-B cell",
  "Dntt+ late pro-B cell"
  )
  
  sub.cell_ontology_class <-
  c(
  "late pro-B cell",
  "late pro-B cell",
  "late pro-B cell",
  "late pro-B cell",
  "late pro-B cell"
  )
  
subtissC = stash_annotations(subtissC, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class)
tiss = stash_subtiss_in_tiss(tiss, subtissC)
data.frame(sub.cluster.ids, sub.cell_ontology_class, sub.free_annotation)

head(tiss@meta.data)

Checking for batch effects

Color by metadata, like plate barcode, to check for batch effects.

TSNEPlot(object = subtissC, do.return = TRUE, group.by = "plate.barcode")

Final coloring

Color by cell ontology class on the original TSNE.

TSNEPlot(object = subtissC, do.label = TRUE, do.return = TRUE, group.by = "free_annotation")

TSNEPlot(object = subtissC, do.label = TRUE, do.return = TRUE, group.by = "cell_ontology_class")

Subset D == hematopoietic precursor cells (Cluster 0)

This subset contains hematopoietic multipotent progenitor (MPP) cells and a very minute fraction of Hematopoietic stem cells (HSC).

subtissD = SubsetData(tiss, ident.use = c(0))

subtissD <- subtissD %>% ScaleData() %>% 
  FindVariableGenes(do.plot = TRUE, x.high.cutoff = Inf, y.cutoff = 0.5) %>%
  RunPCA(do.print = FALSE)

[1] "Scaling data matrix"

  |                                                                                                  
  |                                                                                            |   0%
  |                                                                                                  
  |============================================================================================| 100%

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

PCHeatmap(object = subtissD, pc.use = 1:3, cells.use = 500, do.balanced = TRUE, label.columns = FALSE, num.genes = 8)

PCElbowPlot(subtissD)

sub.n.pcs = 20
sub.res.use = 1
subtissD <- subtissD %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, 
    resolution = sub.res.use, print.output = 0, save.SNN = TRUE) %>%
    RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=65)
TSNEPlot(object = subtissD, do.label = T, pt.size = 1.2, label.size = 4)

previous_annotation_table(subtissD)

[1] "previous_free_annotation"
                                             
                                                0   1   2   3   4   5   6
  common lymphoid progenitor                    0   0   1 123   3   0   0
  granulocyte monocyte progenitor cell          0   0   0   0   2  60   0
  hematopoietic stem cell enriched fraction    41   1   2   0   0   0   1
  HSC enriched fraction                         1  92   0   0   0   0   0
  megakaryocyte-erythroid progenitor cell       0   0   0   0   0   0  46
  NA                                           38   1  21  23  15   4   4
  Slamf1-negative multipotent progenitor cell 188  81 129   9  37   0   2
  Slamf1-positive multipotent progenitor cell   0   0   9   1  78   0   1
[1] "previous_cell_ontology_class"
                                      
                                        0  1  2  3  4  5  6
  granulocyte monocyte progenitor cell  0  0  0  0  2 60  0
  monocyte                              0  0  1  0  0  0  0
  NA                                   38  1 21 23 15  4  4

Check expression of genes of interset.

{r} # subtissD.markers <- FindAllMarkers(object = subtissD, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) #

{r} # subtissD.markers %>% group_by(cluster) %>% top_n(6, avg_logFC) #

The multitude of clusters of each type correspond mostly to individual animals/sexes.

table(FetchData(subtissD, c('mouse.id','ident')) %>% droplevels())

        ident
mouse.id   0   1   2   3   4   5   6
  3_10_M  19 166   5  25  60  45  27
  3_38_F   0   0   0   0   0   1   2
  3_39_F   4   2   2   7   4  13   6
  3_8_M  138   2  83  66  36   4   7
  3_9_M  107   5  73  58  35   1  12

sub.cluster.ids <- c(0, 1, 2, 3, 4, 5, 6)
sub.free_annotation <- c(NA, NA, NA, NA, NA, NA, NA)
sub.cell_ontology_class <- c(
  "Slamf1-negative multipotent progenitor cell",
  "hematopoietic precursor cell",
  "Slamf1-negative multipotent progenitor cell",
  "common lymphoid progenitor",
  "Slamf1-positive multipotent progenitor cell",
  "granulocyte monocyte progenitor cell",
  "megakaryocyte-erythroid progenitor cell"
  )
subtissD = stash_annotations(subtissD, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class)
tiss = stash_subtiss_in_tiss(tiss, subtissD)
head(tiss@meta.data)

data.frame(sub.cluster.ids, sub.cell_ontology_class, sub.free_annotation)

head(subtissD@meta.data)

Checking for batch effects

Color by metadata, like plate barcode, to check for batch effects.

TSNEPlot(object = subtissD, do.return = TRUE, group.by = "plate.barcode")

Final coloring

Color by cell ontology class on the original TSNE.

# TSNEPlot(object = subtissD, do.label = TRUE, do.return = TRUE, group.by = "free_annotation")
TSNEPlot(object = subtissD, do.label = TRUE, do.return = TRUE, group.by = "cell_ontology_class")

Final coloring

Color by cell ontology class on the original TSNE.

TSNEPlot(object = tiss, do.label = TRUE, do.return = TRUE, group.by = "free_annotation")

TSNEPlot(object = tiss, do.label = TRUE, do.return = TRUE, group.by = "cell_ontology_class")

Subcluster the clusters 7 of tiss

Subset E == granulocyte monocyte progenitor cell, monocyte

subtissE = SubsetData(tiss, ident.use = c(7))

subtissE <- subtissE %>% ScaleData() %>% 
  FindVariableGenes(do.plot = TRUE, x.high.cutoff = Inf, y.cutoff = 0.5) %>%
  RunPCA(do.print = FALSE)

[1] "Scaling data matrix"

  |                                                                                                  
  |                                                                                            |   0%
  |                                                                                                  
  |============================================================================================| 100%

Calculating gene means
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating gene variance to mean ratios
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

PCHeatmap(object = subtissE, pc.use = 1:3, cells.use = 170, do.balanced = TRUE, label.columns = FALSE, num.genes = 8)

PCElbowPlot(subtissE)

sub.n.pcs = 20
sub.res.use = 1
subtissE <- subtissE %>% FindClusters(reduction.type = "pca", dims.use = 1:sub.n.pcs, 
    resolution = sub.res.use, print.output = 0, save.SNN = TRUE, force=TRUE) %>%
    RunTSNE(dims.use = 1:sub.n.pcs, seed.use = 10, perplexity=60)
TSNEPlot(object = subtissE, do.label = T, pt.size = 1.2, label.size = 4)

previous_annotation_table(subtissE)

[1] "previous_free_annotation"
    
     0 1 2 3
  NA 3 5 1 2
[1] "previous_cell_ontology_class"
                                      
                                        0  1  2  3
  granulocyte monocyte progenitor cell  0  0  3 66
  monocyte                             99 87 68  2
  NA                                    3  5  1  2

Check expression of genes of interset.

{r} # subtissE.markers <- FindAllMarkers(object = subtissE, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25) #

{r} # subtissE.markers %>% group_by(cluster) %>% top_n(6, avg_logFC) #

The multitude of clusters of each type correspond mostly to individual animals/sexes.

table(FetchData(subtissE, c('mouse.id','ident')) %>% droplevels())

        ident
mouse.id  0  1  2  3
  3_10_M 22  2 30  4
  3_38_F  2 58  6 32
  3_39_F 37  7  4 22
  3_8_M  28  6 20  6
  3_9_M  13 19 12  6

sub.cluster.ids <- c(0, 1, 2, 3)
sub.free_annotation <- c(
  NA, NA, NA, NA
  )
sub.cell_ontology_class <-
  c(
  "monocyte",
  "monocyte",
  "monocyte",
  "granulocyte monocyte progenitor cell"
  )
subtissE = stash_annotations(subtissE, sub.cluster.ids, sub.free_annotation, sub.cell_ontology_class)
tiss = stash_subtiss_in_tiss(tiss, subtissE)

head(tiss@meta.data)

Checking for batch effects

Color by metadata, like plate barcode, to check for batch effects.

TSNEPlot(object = subtissE, do.return = TRUE, group.by = "plate.barcode")

Final coloring

Color by cell ontology class on the original TSNE.

# TSNEPlot(object = subtissE, do.label = TRUE, do.return = TRUE, group.by = "free_annotation")
TSNEPlot(object = subtissE, do.label = TRUE, do.return = TRUE, group.by = "cell_ontology_class")

Final coloring - entire tissue!

Color by cell ontology class on the original TSNE.

TSNEPlot(object = tiss, do.label = TRUE, do.return = TRUE, group.by = "free_annotation")

TSNEPlot(object = tiss, do.label = TRUE, do.return = TRUE, group.by = "cell_ontology_class")

Save the Robject for later

filename = here('00_data_ingest', '04_tissue_robj_generated', 
                     paste0("facs_", tissue_of_interest, "_seurat_tiss.Robj"))
print(filename)

[1] "/Users/olgabot/code/tabula-muris/00_data_ingest/04_tissue_robj_generated/facs_Marrow_seurat_tiss.Robj"

save(tiss, file=filename)

# To reload a saved object
# filename = here('00_data_ingest', '04_tissue_robj_generated', 
#                      paste0("facs", tissue_of_interest, "_seurat_tiss.Robj"))
# load(file=filename)

Export the final metadata

Write the cell ontology and free annotations to CSV.

save_annotation_csv(tiss, tissue_of_interest, "facs")